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Abstract. In view of recent experiments in ultracold atomic systems, the 
photoassociation of Efimov trimers, composed of three identical bosons, is studied 
utilizing the multipole expansion. We study both the normal hierarchy case, 
where one-body current is dominant, and the strong hierarchy case, relevant 
for photoassociation in ultracold atoms, where two-body current is dominant. 
For identical particles in the normal hierarchy case, the leading contribution 
comes from the s-mode operator and from the quadrupole d-mode operator. 
The s-mode reaction is found to be dominant at low temperature, while as the 
temperature increases the d-mode becomes as significant. For the strong hierarchy 
case, the leading contribution comes from a 2-body 5-wave J operator. In both 
cases log periodic oscillations are found in the cross section. For large but finite 
scattering length the amplitude of the oscillations becomes larger in comparison 
to infinite scattering length case. We apply our theory to photoassociation of ^Li 
ultracold atoms and show a good fit to the available experimental results. 

PACS numbers: 31.15.ac,67.85.-d,34.50.-s 
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1. Introduction 

When the properties of a few-body or many-body 
system are insensitive to the details of the microscopic 
interaction between its constituents, the system is 
said to be Universal. A free gas is the simplest and 
somewhat trivial example of universality. A very rich 
and important case is a system governed by the low 
energy scattering parameters. Such scenario can be 
realized in ultra cold atomic systems, where the two- 
body scattering length can be tuned by a Feshbach 
resonance to be much larger than any other length scale 
in the system [HE]. 

The trimer, or the three body system, attracts 
special attention as the simplest non-trivial universal 
system. Moreover, in the 70’s Efimov predicted that 
in the limit of a resonant 2-body interaction, the 
system reveals universal properties mm- A peculiar 
prediction is the existence of a series of giant three 
body molecules, known as Efimov trimers, that was 
verihed experimentally a few years ago. For a review 
see Ref. |S] and references therein. 

The coupling of a radiation field with the system 
is a useful experimental tool in physics. The 
photoassociation and photodissociation processes are 
sensitive ways to measure spectrum and also other 
properties of few and many-body systems [5] . Recently, 
trimers were formed by radio-frequency (rf) excitations 
in both fermionic ®Li [3i and bosonic "^Li |9] systems. 

In a previous work m, we have presented the 
multipole analysis of an rf association process binding 
a molecule of N identical bosons. We have considered 
two different scenarios relevant to photoassociation 
experiments in neutral ultra cold atoms, confined in 
a strong magnetic held: (i) spin-hip process, and (ii) 
frozen-spin process. 

In the spin-hip mechanism, the photon induces a 
spin state change in the absorbing atom, and most of 
its energy is dedicated transferring the atom from one 
hyperhne state to another. The typical energy of the 
photon in this case is of the order of the hyperhne 
splitting. The spin-hip mechanism was materialized 
experimentally by m iHi- Previous analysis of rf 
experiments dH Ha [I3II1 US], which relied on the 
Franck-Condon factor, is appropriate for describing 
these spin-hip reactions. 

In the second, frozen-spin, mechanism, the photon 
energy is insufficient to bridge the hyperhne gap and 
therefore the atom spin is “frozen” throughout the 


2 

process. In this case the photon induces a molecular 
bound-free transition and the typical photon energy is 
of the order of the molecular binding energy. 

In |10j we have shown that the spin-hip and 
frozen-spin processes differ by their operator structure 
and by the de-excitation modes that contribute to 
the photoassociation rate. The photoassociation 
experiments carried out by the Bar-Ilan group, see e.g. 
|S], relied on the frozen-spin mechanism. 

We have studied one-body current in frozen-spin 
reactions, dealing with dimer formation |10j . and also 
calculated numerically the quadrupole response of a 
bound bosonic trimer m- Recently we have identihed 
another scenario, where two-body currents dominate 
the frozen-spin process, and explored the occurrence 
of log-periodic oscillations in the reaction cross section 

m- 

Here, we study few more aspects of the three body 
photoassociation process. These are (i) The interplay 
between the s-mode and the d-mode reaction channels 
at different temperatures and photon energies, (ii) 
The effect of finite scattering length on the reaction 
rate. We utilize the hyperspherical coordinates 
with the adiabatic expansion to study the trimer 
photoassociation process, for both one-body process 
as well as for the two-body process. Analytic results 
for the transition rates at the unitary point (where 
the scattering length diverges) are derived using the 
zero-range approximation, and numerical calculations 
complete the picture for a large but finite scattering 
length. Similarly to the dimer case [TUI, the s- 
mode and the d-mode are found to be the leading 
order contributions in the one-body process. The 
relative importance of these modes depends on the 
temperature, where at low temperature the s-mode 
is dominant while at higher temperature the d-mode 
becomes as significant. For both one-body and two- 
body processes, log-periodic oscillations modulate the 
cross section. We found that comparing with the 
unitary point, the amplitude of these oscillations is 
somewhat larger for finite scattering lengths. 

As a concrete example for the usefulness of 
our approach we analyze the dimer and trimer 
photoassociation experiments carried out by the Bar- 
Ilan group in an ultracold atomic ^Li gas [Uj. As 
already pointed out, the photoassociation mechanism 
in these experiments is the frozen spin process. Due 
to the strong hierarchy, see Sec. the coupling of 
the photon to the system is dominated by the 2-body 
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magnetization current. Our theory fits well to the 
available experimental results. 

The paper is organized as follows, in section 2 
we introduce our model and reiterate the multipole 
and the currents expansions, describing the interaction 
of an rf field with the atomic system. In section 

3 we introduce the three-body problem and solve it 
using the hyperspherical adiabatic expansion. The 
transition matrix elements are calculated in section 

4 for infinite, and also for large but finite, scattering 
lengths. The transition rates are calculated in section 
5, and a detailed comparison to experiment is presented 
in section 6. Section 7 concludes our study. 


2. The model 


In photoassociation experiments [71 m IS] , a few MHz 
rf radiation field is applied to the ultracold atomic gas. 
Stimulated emission occurs when the photon energy 
matches the difference between the ultracold gas 
atoms and the molecule energy, resulting in molecule 
formation. The transition rate of such process is given 
by Fermi’s golden rule, 

-Ef- hcv), (1) 

* / 

where an average on the appropriate initial 

continuum states and is a sum on the final bound 
states. The coupling between the neutral atoms and 
the radiation field takes the form 

Hi = — dx/i(x) • V X A(x) 


where /i is the magnetization density and A is the 
electromagnetic field at position x. We have shown 
that in effective low energy theory m H has 
contributions from one-body current as well as more 
body currents 

M(x)=/r«(x) + /r(2)(x) + .... (2) 


An effective theory has some UV cutoff A, reflecting 
some short-range physics which is ignored. Naive 
power-counting suggests that each order in this 
expansion is suppressed by a factor of (Q/A)^, where 
Q is the typical momentum of a particle in the system 
under consideration. The leading order (LO) one-body 
current is m 




(1) 


(x) = ^Sj/ri(5(x-rj), 


(3) 


where fXi is the magnetic moment of a single particle, Vj 
and Sj are the position and the spin of particle j. The 
two-body contribution to the atomic electro-magnetic 
current enters at the next order (N^LO), or (Q/A)^, 
and takes the form [19] 




i<j 


-)(5A(r,-r^) (4) 


where L 2 = L 2 W is the coupling constant between 
the radiation field and the four boson fields. The 
notation i5A(r) stands for Dirac’s 5-function smeared 
over distance Ti/A. 

In the low energy limit Q <C A, more particles 
current can be ignored due to suppression by additional 
(Q/A)^ factor associated with each extra particle field. 

Using box normalization of volume H, the electro¬ 
magnetic field reads 



where C = 1,2 are the linear photon polarizations, w is 
the photon frequency, k its momentum, and Cq is the 
vacuum permeability. 

For bosonic systems confined in a strong (static) 
magnetic field, the initial and final state atomic wave 
functions can be written as a product of spin and 
configuration space terms, d' = iPlmXMip- The 
spin wave function XMp is the symmetrized function 
XMf = S [|tof(1))|top’( 2)) ... |m_F(A))] with magnetic 
quantum number Mp = mp^j). The spatial wave 
function tppM = V’iM(ri, r 2 ,... tat) is a symmetric 
function with angular momentum quantum numbers 
EM. Using this factorization, the one-body transition 
matrix element in Q takes the form 

(/.KIhUi.) = 

N (5) 

1=1 

whereas the two-body part reads 


(/,kC|i7r^|z) 


• I ^ ^2 

V 2Xtujeo A^ 


N 

|(s, -f Sj) • (k X ekc)IXMp)^ 


i<j 


( 6 ) 


Using spherical notation, the spin operator reads 
s-(kxekc) = EA(“)^'®-A'(kxekc)A, where A = 0,±1. 
The geometrical factor, rj = (k x ekf)A, is maximal 
(minimal) for the frozen spin A = 0 case (spin flip 
A = ±1 case), when the rf magnetic component is 
parallel to the static magnetic field. 

The photon wavelength of rf radiation is much 
larger than the typical dimension of the system R, 
therefore kR ^ 1 and the lowest order in kR 

dominates the interaction. When the photon can 
induce a Zeeman state change, i.e. spin-flip, the leading 
contribution comes at order k. Energy is delivered to 


2 
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the system through the spin matrix element and we 
can approximate ~ 1, to get 

K/.k<iHri.)p=^ 


2r2ceo 


E 

A=±l 


N 


(7) 


I (Xm' I E n IV'Im) ?h,L'5M,M' ■ 




The last term on the rhs of Q is just the Franck- 
Condon factor, that appears in the commonly used 
theory of photoassociation, such as Eq. (3) in Ref. 

m- 

When the rf photon cannot induce change in the 
spin structure of the system, A = 0, and we call the 
process a frozen spin reaction. In this case Mp = Mp, 
IXmjt) = IXmjt) fh® transition matrix element can 
be written as 


N 

i=i 

where (so) = ^ Ej (Xmj. Isj.oIxmj.) is th® average 
single particle magnetic moment, which plays the role 
of an effective charge. 

In the long wavelength limit, the exponent can be 
expanded to yield 


N 

E' 


N 


,2k-r^ 


N 


N + tJ2k-r,--y,k 


2^2 

3 






- ^ E E Y2-mik)Y2m{fj), 


N 


15 


j=l m 


where Yim are the spherical harmonics. Clearly this 
expansion has transparent physical meaning. The zero 
order operator is proportional to 1 and stands for 
elastic interaction. In photon emission reaction this 
process is forbidden by energy conservation. Next 
comes at first order the dipole. Dealing with identical 
particles, this term is proportional to the center of mass 
and hence cannot affect internal degrees of freedom. 
Two operators appear at second order: the operator, 
corresponding to s-mode reaction, and the quadrupole 
terms, corresponding to d-mode reaction. Summing 
over the initial and final magnetic numbers M and M', 
the transition matrix element reads 1101 


E 


(i&)i 


M,M' 


A'Ktikf’iJL{rf 
2Dcen 


i(so)r 


N 


i=i 1=1 


For the two body current, we can again utilize the 
long wave approximation r:; l to get 

Yhk L\rf 


i(/,kciidr)i*)p= 


i(so)r 


KV’l 


2r2ceo A® 

( 11 ) 

i<j 

To find the most relevant operators for the 
photoassociation process one should consider both the 
low energy expansion (|^ and the long wavelength 
expansion Dealing with trimer photoassociation, 
the photon energy has to be of the order of the 
trimer binding energy E 3 . Therefore the long 
wavelength expansion parameter can be written as 


kr 


yj E^/Mc 


Q/Mc, where we used r 


Ti^/ME^. Therefore one should compare the long 
wavelength expansion parameter, kr « Q/Mc to the 
low energy expansion parameter Q/K. If Q/A ^ 
K/Mc, the two-body currents appearing at order 
(8) (Q/A)3 are much smaller than the second order 

(A;r)2 « {Q/Mcf^ terms. This normal hierarchy case 
is the situation in the limit Q — > 0. In the other 
extreme, when A Me, the two-body current is more 
important than the one-body current, proportional to 
(Q/Mc)^. This case of strong hierarchy is typical 
for photoassociation experiments in ultracold atoms 
[3 mg. There, the separation between the mass scale 
and the binding energy is much bigger than the ratio 
between the scattering length and the effective range. 


(9) 3. The Three Body Problem 


The Schroedinger equation governs the dynamics of a 
quantum 3 particle system 

{T + W)f3 = Ei}, (12) 

where T is the kinetic energy operator in the center of 
mass frame and W is the potential. Here we assume 
a short range 2-body forces, thus W = J2i<j ^(ki ~ 
rjl). To solve this problem we use the hyperspherical 
coordinates and the adiabatic expansion [21111]. In 
the following section we review these methods and 
extend them for the L > 0 case. 

3.1. The Hyperspherical coordinates 

To eliminate the center of mass motion, we define the 
Jacobi coordinates, 

' rk 


= \ 0(^3- T^k), 


y* = \/ 3 1 -D + 


2 v 3 V 3 V “ ' 2 

vhe^-e^^{ijfc} is a cyclic permutation of {123}. The 
'^r^lation between two sets of Jacobi coordinates (xi,yi) 
id (Xj, Yj) is given by the kinematic rotation. 


— cos (pij sin (j 3 i 

— sintii,- —cosd 


Xi 

yi 


( 13 ) 
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For the case of 3 identical particles, (fij = ± 7 r/ 3 , where 
the sign is fixed by the parity of the permutation {ijk}. 
For each set of Jacobi coordinates we can define a set of 
hyperspherical coordinates (p, fli), where = xf +yf, 
and tanoi = Xifyi. Two different sets 
of hyperspherical coordinates are related through 

2 sin^ aj = 1 — cos cos 2ai — 7 ^ sin 20^ sin 2ai, (14) 

where 7 i = Xi ■ yi. The hyperradius p is invariant 
under kinematic rotations, particle permutations, and 
therefore under change of coordinates sets. 

In the hyperspherical coordinates the kinetic 
energy operator reads. 


T = - 


2m i 9p2 p dp 


(15) 


Here, is the square of the grand angular momentum 
operator, 


= - 


1 


92 


sin 2 q! da^ 


sin 2 a 


ll 


sin^ a 


+ 


‘■y 


-4, (16) 


and lx, ly are the angular momentum operators 
corresponding to x, y coordinates. The hyperspherical 
presentation of a central two-body potential reads, 


V{\r^ - Tjl) = y^V{V2psm aj). 


i<j 


3.2. The Adiabatic Expansion 


(17) 


Next we apply the adiabatic expansion | 20 j and write 
the wave function in the form 


^(p,H)=5]p-5/27^„(p)$„(p,H). 


(18) 


The hyperspherical functions $„(p, H) are the solutions 
of the hyperangular equation, 

(^K'^{V2p sin a,)+ ^ ^n{p,^) = 

r'2$„(p,H), 

corresponding to the eigenvalue u^. The hyperradial 
functions TZn{p) are the solutions of the hyperradial 
equation, 

+ ^eff(p) - 7^„(p) = 

{‘2Pnn'-p, - \- Qnn')'P-n'{p) 

,,, 9p 

n';n'^n 

where e = 2mE/h^, t4ff is the effective hyperradial 
potential, and Pnn', Qnn' are the non-adiabatic 
couplings. The effective potential is given by 


( 20 ) 


Ves[p) = 


^lip) - 1/4 


Qnn, 


and the non-adiabatic couplings are 
9 


Pnn'(.p) — \ 4*n(P;fl) 


Qnn'ip^ — \ 4in(P;fl) 


dp 

92 


dp"^ 


(P; fl) 


4>n'(p, fl) 


( 21 ) 


The expectation value (...)n stands for integration over 
the hyperangles H. 

3.3. The Hyperangular Equation 

For low energy physics, when the extension of the 
wave function is much larger than the range of the 
potential, one can utilize the zero range approximation. 
In this approximation the lateral extension of the 
potential is neglected all together, and the action of 
the potential is represented through the appropriate 
boundary conditions. For a two-particle system the 
low energy interaction is dominated by the s-wave 
scattering length a and the wave function fulfills 
the boundary condition [u'/u]r=o = —Ija. The 
corresponding 3-body condition is 

^ ^ ^ =-V2P. (22) 

a 


2 a, 4) da. 


-2a,^ 


cti—0 


Dealing with wave functions of definite total 
angular momentum quantum numbers L,M, we use 
a Faddeev-like decomposition, 

4>(p,H) =N,Y.Y. (23) 

^ lx 

where Ni is a normalization constant and 

mx,my 

It should be noted that for a bosonic system U must 
be even. 

In the zero range approximation each Faddeev 
component is a solution of the free hyperangular 
equation, 

{k^ - A)cf,{p,^,)YY^{x,y) = 0 , (25) 

with the appropriate boundary conditions. 


One family of solutions to (25) is the hyperspher¬ 
ical harmonics 1221. 


where 

pf i i sin^"’acos*“ aP,^*‘''''^/2J«+i/2)('(.Qg2Q,).(27) 


(26) 


Here P„ 


1/3) 


are the Jacobi polynomials, n = {K — lx — 


ly)/2, and 


Nki = 

lx Hy 


I 2ti\{K + 2)(n lx 1)! 

r(7i lx 3/2)r(?i ly 3/2) 
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Note that n is a non-negative integer and therefore 
K > U + ly. The hyperspherical harmonics correspond 
to A = K{K+A) and are regular at a = 0 and a = 7r/2. 
Another family of solutions can be found by 


generalizing (26) to a non-integer K (and n ), 


loosing the regularization in one edge. Due to the 


boundary condition (221, we need such solution for the 


interaction channel, i.e. the U — 0 component. Such 
solutions, corresponding to the eigenvalue K = — A, 

are P^{a )/sin 2a and Q^{a)l sin 2a [23], where 


P^ia) = cos^ a 
and 


/ d 1 

(a) = cos'^ a 


(^ 



1 5a 

- Sin 1 

/ ( a — — 

cos a J V 

V 2// 


sin (va). 


(28) 


(29) 


where the last term is the 3-j symbol. Some conclusions 
emerge: 

a. The projection at a^ = 0 is zero for all waves 
but U = 0,ly = L; similar considerations show that 
the projection at a^ = 7r/2 is zero for all waves but 

lx — T, ly = 0. 

b. Due to the (t)^=' factor, for odd the two 
rotations cancel each other. 

c. The 3-j symbols equal zero unless I'x + I’yL is 
even. 

d. Using this result, an explicit formula for the 
lx = Q or ly = Q Raynal-Revai coefficients {lxly\l'xl’y )kl 


[24] can be constructed. See Appendix A for details. 


\ da cos a 

The first (second) solution is regular at a = 7r/2 
(a = 0) and not regular at a = 0 (a = 7r/2). 

For evaluating the rf transition matrix elements 
we need to consider L = 0 and L = 2 states. Their 
explicit form is given by 

P»=si„(„(a-|)) 

Q°(a) = sin(j^a), 
and 

P^{a) = 3^ tan a cos ^a — 

— {2 + v‘^ — 3 sec^ a) sin ^a — ^ 

Q^{a) = 3^tanacos(z^a) — {2 + — 3sec^ a) sin(j^a). 

3.4- The Projection Operator 


Focusing on the lx = 0, ly = L case, the projection 
operator can be computed for any a^ [2S]. This 
is because for a^ < 7r/3 the rotated function is an 
eigenfunction of the kinetic energy operator with the 
same eigenvalue, but is regular for a^ = 0, therefore 
it is proportional to Q{j. The proportionality constant 
can be calculated by matching with the known value 
at ai = 0. Similarly, the rotation for a^ > tt/3 is 
proportional to P^, and the proportionality constant 
can be found by matching at a^ = 7r/3. 

Finally, defining P^ = P^/sin(2a) and = 
Q^/sin(2a), 


P 


(0L)(0L) r pL 


[Pj^(a,)](a.) = (-l)^x 


(31 


Q^(a,)Pi'(V3)/Q^(0) 

P^(az)Q^(^/3)/g^(0) 


0 < ai < 7r/3 (35) 

7r/3 < ai < 7r/2. 


It is convenient to project the wave function (23) on a 
single coordinate system |23j . The resulting expression 
takes the form 




3.5. Imposing the Boundary Condition 

In the limit of infinite scattering length, the adiabatic 
expansion decouples and the non-diagonal couplings 
Pnn'TQnn' Vanish [23] . In this limit, the lowest energy 
hyperangular function takes the form 

(jiiai) = P^(ai)/sin(2ai) . 


^xi^y 

Y. E («.)}. 

where Rij, the projection operator is given by 
dPyWj'y) 

(33) 


(36) 

^o^The eigenvalues v will emerge as we impose the 
appropriate boundary conditions at ai = 0. We have 


seen in (34), that the rotation Rij at a^ = 0 includes 


R. 


[(j){aj)\{ai) = 

dx, J dyiYi^^*{x„yi)(j){aj)Yifi¥{xj,yj) 


only the lx = 0 partial wave, therefore for small ai (32) 
is simply given by [26] . 

=Pi;{a.) + 2(1 + (-)^)a.P,^(7r/3) 

YO{ai) 


To calculate R^ we first study the limit a^ = 0. 
From (13) it follows that in this case, Xj = ±-\/3yi/2. 


Yj = —Yi!^ and aj = 7r/3. The ± sign depends on 
the direction of the kinematic rotation. The integral 


Therefore the boundary conditions 
dP^ifS) 4 


turn into 


da 


^(l + (-)^)Pj^(V3) = -^PiifS) .(38) 


in (33) can now be evaluated, yielding 


For L = 0 the resulting equation for z/ is [3] 


R^^umy)^^^aMo^, = 0) = 

cj{7r/3){TfSi^,oSi,,LyJ{2l'x + 1)(2Z; + 1) 


8 


I'. ly L 
0 0 0 


I/CQf 

(341 


V2p . 


Vcos{v'K/2) -— sin(i/7r/6) = - - sin(zz7r/2). (39) 


73 
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Figure 1. (Color online) The two lowest angular eigenvalue 
as a function of p/a for L = 0 (Upper panel) and L = 2 
(Lower panel). The asymptotic solutions are dashed. Negative 
1 / corresponding to imaginary solution. 


For |a| = oo the solution with lowest is vq « 

1.00624/, corresponding to the Efimov trimer. For a > 
0 this solution approaches asymptotically a particle 
scattering from a universal dimer, where Ves{p — 
00 ) = — l/o^ [27]. The spurious solution u = 4 is just 
$ = 0 . 

For L = 2 the corresponding equation for v is, 
i/{A — v^) cos{v'k/2) + 24i/cos(;/7r/6)+ 

— 10) sin(z^7r/6) = —— 1) sin(z/7r/2) ^ ^ 

V 3 o 

For |a| = oo the lowest non-trivial solution is V 2 ~ 
2.82334. For a > 0 this solution asymptotically 
converges to particle-dimer d-wave scattering, where 
Ves(p —> oo) = L(L + l)/p^ — 1/a^. The ^ = 0,1, 2 
correspond to the spurious solution $ = 0. For p |a| 
solutions with v an even integer, are just the regular, 
free, hyperspherical harmonics. 

The two lowest eigenvalues of the hyperangular 
equation, with the appropriate boundary conditions 
), are plotted in Fig. [^for L = 0 and L = 2, 
with their asymptotic forms. 


(391,(40 


3.6. The Hyperradial Equation 


To proceed analytically, we focus on the unitary limit, 
|a| —>■ oo. In this case vl{p) = vl and <l>(p, fl) = $(11), 
therefore the non-adiabatic couplings, Eq. (21), vanish 


and the hyperradial equation for TZ{p)/^ is just the 
Bessel equation, 

d^TZnip) , - 1/4, 


dp^ 




-TZn{p) = UP-n{p)- 


(41) 


In order to calculate the trimer photoassociation, 
caused by s— and d-wave operators, we have to 
solve the L = 0 bound state as well as the L = 
0, 2 continuum states. The needed solutions are the 
followings: 

I. A bound state (e < 0) with L = 0, and 
Vq « 1.00624/. In this case the relevant solution is 
proportional to the modified Bessel function of the 
second kind and imaginary order, y/pK^g^np), where 
K = \f-^. At the origin, this solution behaves like 
sm.{v\n{Hip/2) — where « —0.301. Therefore 
regularization is needed to avoid collapse, e.g. setting 
Ti-{p < Po) = 0 for some finite po- The result is the 
discrete Efimov spectrum. 


515“ 


^ g-277m/|iyol 

Co 

The normalized wave functions are 


^^\p) = NBKmy/pKi,g{KmP) 


where Nb = \/2 sin vq-kIvq-k « 2.730, and 
22.7“"". 


(42) 


(43) 


(44) 


KO 

II. Scattering state (e > 0) with L = 0, and 
r'o « 1.00624/. The solution is composed of the real 
part of the Bessel functions of the first and second kind 
of imaginary order, 

T^si.p) = [sin5Re[J^o(gp)] -f cos(5Re[l/.o(gp)]] (45) 


where q = y/e and we assume normalization in a sphere 
of radius R. The phase shift <5 is to be found from the 
boundary condition, TZs{pq) = 0. 

In the scattering problem, functions with higher 
v, corresponding to the same angular momentum 
L = 0 and energy e, are also legitimate solutions. 
However, due to the orthogonality in fl, there is no 
overlap between these functions and the bound state, 
and therefore no contribution to the photoassociation 
transition matrix elements. 

III. Scattering state (e > 0) with L = 2, and real 
1 ^ 2 . The solutions are composed of the Bessel functions 
of first and second kind, 


'^d(p) = Y ^ [sin SdJi^g (qp) cos SdY^g {qp)] . (46) 

The phase shift 5d is to be determined by the condition 
T^diPo) = 0- The three lowest values for 1^2 are, 
1^2 « 2.823, « 5.508 and « 6.449. 
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4. The Transition Matrix Elements 

Now that we have solved the Schroedinger equation 
and obtained the L = 0 bound state wave function 
and the L = 0,2 scattering wave functions, we are in 
a position to evaluate the transition matrix elements, 
in both hierarchies. For the normal hierarchy case the 


needed matrix elements (101 read. 


Kb = 


N 

^ ' 3 
i=i 


of Isin, q) is obtained when q « k/2. For higher Efimov 
states, this integral scales such as 

p^^'^Isium.qp) = IsiKo,q) 

where p = Kml^o = 22.7“™, 

It should be noted that (521 leads to an oscillatory 
log-periodic response function, see discussion in 
m- At threshold, the matrix element ([5^ 




for the s—mode transition and 

i=i 


(47) 


(48) 


gets a particularly simple form which can be well 
approximated by m 

Bs 


Is{K,q) « Ck -^y/q 1 -h 


cos(2|i/o|ln-) 
2 K 


(54) 


for the d—mode. For the strong hierarchy case 0 the 
matrix element is. 


N 


hb = - Yj)\\f}l), 

i<j 

These matrix elements are proportional to integrals of 
the type 


(49) 


(50) 


pOO 

Ij{u,m)= / dpK,y^{Kp)p'^^^Re[J^{qp)]. 

■3 PO 

or Iy(i/, to) with Yi, replacing J,^. 

Taking the lower limit to zero, these integrals can 
be evaluated analytically, 

'2™ AT™ 1 

"1^ (-) 2 Fi{a,b;c-,-{q/Kf) (51) 


where C is a constant that contains the normalization 
factors, and B 3 « 8.475% is the normalized amplitude 
of these oscillations. These oscillations modulate the 
matrix element all the way to the high energy tail. 
Whereas log periodic oscillations at the high energy 
tail appear in all partial waves [28], at threshold the 
oscillations appear only in this s-wave transition. 

The quadrupole matrix element. The quadrupole 
operator, '^jr'^Y^{fj) connects the L = 0 bound 
state with L = 2 scattering states. Using = 

R-cm — 


Yj and working in center of mass coordinate 
system, this operator reads | y?^ 2 ^(di)- Setting 


= p'^cos'^aj, 


3 S] "2 yyj) 
the reduced matrix element can be 


Pj{i 3 , to) « Re 


--m+2 


written as 




3 


where 2^1 is the hypergeometric function with 
parameters a = + 1, 6 = 1, and 

c = V + \, and = F(a)r(6)/F(c). For q < k this 
approximation amounts to an error tYLjX ~ 
where xq « 0.06 is the largest zero of K,^g{x) and Xn is 
the n’th root. Similar results can be obtained for Py- 

4 . 1 . Transition matrix elements in the normal 
hierarchy case 

The matrix element. The operator connects the 
L = 0 bound state to the L = 0 scattering state. 
We note that r'j = p^ + 3 Rqm where Rcm is the 
center of mass coordinate. As the CM cannot induce 
inelastic transition, the matrix element is reduced into 
the hyperradial integral 


Here the hyperradial integral Ip is 

poo 

dpnUp)p^Rdip)- 




(55) 


(56) 


Substituting Eqs. ( |43t and ( |4^ into ( |^ , we get an 
expression similar to (53), 


/ '^Q 

Ip{K,q) = \ —NBK[Ij{iy 2 ,‘ 2 )smSd + lY{i^ 2 ,‘ 2 )cosSd] .(57) 


R 


The hyperangular integral Iq, reads 
In = y dU$/(U) y^cos^ ajY2o{yj)^i{Q). 


(58) 


To evaluate this integral, the hyperangular wave 
function, written in mixed coordinate systems. 


$l(U) = ) ^P,f{ai)YQL{Xi,yi)/sm2ai 


Is{K,q)= dpn*B{p)p^Tls{p). 

Jo 


(52) 


Substituting Eqs. (43) and (45) into (52), we get 


is to be integrated in a single coordinate system. 
Transforming to body-hxed coordinate system and 
integrating over the Euler angles, the integral over the 
R{K,,q) = J ^NBK.[Ij(va,2)sm.5 + Iy{vq,2)cosS\ angles {ai,xi,yi) is transformed into an integral 

V over (ai, 7 i) [30], that is evaluated numerically. See 

This integral depends on the momentum of the bound [Appendix B| for details. For the lowest 1^2 the integral 
state K and the scattering state q. The maximal value yields Iq « 0.368. For the next two the integrals 


q 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 

q/A- 


Assuming a simple regularization, (5 a( x) = 3A^0(1 — 
Aa;)/47r, this matrix element can be calculated, 

hb sin(i/o7r/2)px 


647r 


R' 


(60) 


[I,j{vo, -1) sin (5 + -1) cos (5]. 


To renormalize the transition matrix element and 
remove the cutoff dependence we turn to the dimer 
photoassociation m- In the zero range approximation 
the dimer s-wave continuum wave function is just tpi = 
•\/2/i?sin(gr+(5)/r, and the corresponding bound state 
wave function is 'i/'s = •\/2Ke“”’’/r. Consequently the 
transition matrix element due to the 2-body process is 


Figure 2. (Color online) One-body current: The transition 
matrix element 7^^, as a function of the relative momentum 
normalized to 1 for the maximum value at |a| = oo The two 
modes are presented - (red, peaked around q = and 

quadrupole (blue, peaked around q = k.), with their sum (black). 
Shown are the results for the unitary point (solid lines) and for 
finite scattering length, K.a = —100 (dotted lines) and na = —50 
(dashed lines). 


are respectively 0.023 and 6.9 x 10“^, therefore these 
modes give negligible contribution to the reaction 
rates. 

Here as in the previous case, the hyperradial 
integral Ip leads to log-periodic oscillations in the high 
energy tail m- 

Relative importance. The relative contribution of 
the s, d modes to the trimer formation is displayed in 
Fig. J where the last term in parenthesis on the rhs 
of (101 is presented normalized, along with the s and 
d components. Similarly to the dimer formation case 
HD], the s-wave association is peaked around q = k/ 2, 
while the d-wave association is peaked around q = k. 

In case of large but finite scattering length, the 
eigenvalue of the hyperangular equation, (391 and (401, 
depends on p and therefore the hyperradial equation 
(201 is to be solved numerically. In addition, the 
non-adiabatic couplings (21) must be considered. The 
numerical results for finite a are also shown in Fig. 
for Ka = —100 and na = —50. From the figure it can 
be seen that the |a| —> oo limit is well reproduced if 
kIoI > 100. 


4-2. Transition matrix elements in the strong 
hierarchy case 

Moving now to study the two-body current, there is 
no center of mass contribution, and the transition 
operator is a scalar exciting an L = 0 bound state 
into an L = 0 scattering state. The transition matrix 
element is given by, 

hb = = ^(?/’B|<5A(x)|i/'i) -(59) 


hb,d =(V'B|(5A(i')|V’i) = D 


- 


q cos(d 



q cos S + K sin S 


q^ + 

«;sin((5 -I- g/A)\ 


(61) 




V 


which for large values of A can be approximated as 


hb,d = A^ —W — sin(5 -h 0(A) 


(62) 


Comparing this result to ( |60| ), one can see that the 
cutoff dependence of the 2-body current operator, L 2 = 
L 2 (A) in can be fixed by dimer photoassociation 
experiments such that the trimer photoassociation rate 
is a prediction of this theory. 

As an alternative approach to fix L 2 (A), we can 
study the shift of the dimer binding energy due to a 
small change in the static magnetic field. On one hand, 

SE 2 = (V'sIF^Vb) = ^ —(so)A2(5H + 0(A). (63) 
A'^ ira 

On the other hand, in the vicinity of a Feshbach 
resonance 

a{B) » Ubg ^ ^ , (64) 

where Ubg is the background scattering length, A is the 
resonance width and Bq its position, therefore 
-2 


SE,= ^fbiB = 


2 h^ 


da dB maabg2^ 


SB. 


(65) 


Comparing (63) and (65), T 2 (A) can be determined. 


A2(A) = 


27r 


-A. 


( 66 ) 


3 mobgA (so) 

In Fig. (|^ the two-body transition matrix element 
is shown as a function of the relative momentum q/K. 
The validity of the approximation deriving (51) is 
also checked by evaluating the integral numerically, 
starting from pQ. It can be seen that while for the 
low momentum regime this approximation is excellent, 
for the high momentum tail it deviates from the exact 
result. Numerical results for finite a are also shown in 
Fig. For na = —100 and na = —50, we see no 
significant difference comparing to the unitary point. 
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q/K 

Figure 3. (Color online) Two-body current: The transition 
matrix element as a function of the relative momentum q/k., 
normalized to 1 for the maximum value at |a.| = oo. For the 
unitary point the exact result is shown (blue solid line), as well 
as the approximation of | |60[ l (red dot-dashed line). Also shown 
are results for finite scattering length, Aca = —100 (dotted lines) 
and Ka = —50 (dashed lines). 


5. Transition rates 

To evaluate the trimer photoassociation rate 0, we 
have to average on the initial states and to sum 
over the final states Yf- The initial three-body 
state \i) = \qLMlocly) describes three atoms in the 
continuum with relative momentum q and angular 
momenta L, M, lx,ly The average on these states takes 
the form 

where P[K\{,<i) is the probability of finding an atomic 
trio with quantum numbers [K] = {KLMlxly} and 
momentum q. For system in thermal equilibrium 
with temperature T higher than the condensation 
temperature, 


1 __«V_ 

P[K]iq) = ^e 

where Z is the partition function, ks is the Boltzmann 
constant and T is the temperature. 

The sum Yf contains all possible trimer quantum 
numbers and all possible photons weighted by the 
emission function, 

+ ( 68 ) 

/ k,C L’,M' 

where is the number of photons with momentum 
k and polarization ^ in the final state. The stimulating 
rf radiation is a narrow distribution centered at some 
krf. Therefore iVk^ « Vk,f(;,.,(5k.kjf<5^.Crt j where kj-f is 
the momentum of the stimulating rf field and is its 
polarization. 


5.1. Transition rate in the normal hierarchy case 


Substituting (101, 
rule 0, the trimer formation rate in 
hierarchy case is given by 


( [67| ) and ( |68| into Fermi’s golden 
the normal 


16 inulRm . 


fv^jPiq) 


IJ 

02 


5^ 

152 


.(69) 


where Uri = /fl is the photon density of the 

rf field. The relative momentum q is connected to the 
photon energy through energy conservation h^q^/2m = 
—E^ + hui, where E 3 is the trimer binding energy. 

The relative importance of the s, d modes shifts 
with temperature. This point is demonstrated in Fig. 
1^ where the s, d rates to the trimer photoassociation, 
at the unitary point |a| = 00 and for large but 
finite scattering length are presented for = 

E 3 and ksT = O. 2 E 3 . The importance of the 
d mode grows with the ratio ksT/E^. Therefore 
we can conclude that for small ksT/E^ values, the 
photoassociation is an s-wave process while for large 
values the two modes have similar importance. In 
addition, oscillations in the response can be clearly seen 
at the low frequency regime. These oscillations result 
form the log-periodic structure of the matrix element 
m- At higher temperature these oscillations are 
stronger while at lower temperatures this manifestation 
of the Efimov effect is suppressed by the Boltzmann 
factor. The amplitude of the oscillations is larger for 
finite scattering length than at the unitary point. 


5.2. Transition rate in the strong hierarchy case 


For the strong hierarchy case, the trimer formation rate 
is given by 


2 b,t 


Li, AmR 


2 \(^o)\^d'^-Pt{q)l2b,f 

A"" n Co Q 


OJ 


(70) 


Here we have added the subscript “t” to P{q),l 2 b 
to stress these are trimer related quantities. The 
corresponding dimer formation rate for this case is 
given by 


2 b,d 


L 2 2mR 
V n^c'^eo 


nrf \{so)\^ ri^-Pd{q)llb,d^ 


(71) 


where Pd{q) is the probability for a pair of particles to 
be in L = 0 state with momentum q. Comparing these 
two results and using Eqs. (62),(66l we finally obtain 

,Pt{q) ^ 2 b,t 


2 b,t 


2 b4 


= 2 


Pd{q) III 


(72) 


26 ,d 


and 


2 b,d 


mobgA 


2 


2 m 

Lc^cq 


nrfv'^ 


-Pd{q)KsiiL 6 (q) (73) 
9 


One should note that these results are independent 
of the cutoff parameter A. The dependence of the 

















Photoassociation of Universal Efimov Primers 


11 




0.001 0.01 0.1 1 10 





Figure 4. (Color online) One-body current: The normalized 
trimer photoassociation rate as a function of the rf photon 
frequency, for different temperatures and scattering lengths. The 
total (black), the (red, low energy peak) and the quadrupole 
(blue) rates are presented, for temperature of ksT = Ez (upper 
panel) and k^T = 0.2Ez (lower panel). Shown are results for 
the unitary point (solid lines) and for finite scattering length, 
Ka = —100 (dotted lines) and na = —50 (dashed lines). 


rates on the normalization radius R enters through the 
probability Pd/t{q)- 

The trimer transition rate for the strong hierarchy 
case is shown for ksT = in Fig. Here as 
in the normal hierarchy case, log periodic oscillations 
are visible. For finite scattering length they are 
amplified. Since there is only a single channel here, 
the low temperature behavior is similar to the high 
temperature, and therefore not shown. 

6. Experimental results 

In this section we apply our theory to the photoasso¬ 
ciation of ultracold ^Li atoms, and fit it to the experi¬ 
mental results of the Bar-Ilan group [5]. 

In this experiment, a gas of ^Li atoms on the |F = 
l,mF = 0) state was cooled down to T « with 


Figure 5. (Color online) Two-body current: The normalized 
trimer photoassociation rate as a function of the rf photon 
frequency. Shown are the results for k^T = Ez at the unitary 
point (solid lines) and at finite scattering lengths, na = —100 
(dotted line) and K.a = —50 (dashed line). 


mean density of about 10^^cm“^. A static magnetic 
field was applied in the vicinity of a Feshbach resonance 
at ^894 G. Then an rf field was turned on for ~ 1 msec. 
The magnetic component of the rf field was parallel to 
the static magnetic field, therefore ry = 1. The number 
of atoms remains in the trap was measured as function 
of the rf frequency. 

The scattering length in this experiment is 
positive, therefore both dimers and trimers can be 
associated by the rf field. The experimental signal 
for photoassociation is particle loss from the trap. In 
both cases three particles are lost for each associated 
molecule. Internal collisions break the trimer into a 
deeply bound dimer and an atom, both carry high 
kinetic energy and therefore escape the trap. The 
dimer collides with a third atom and decays to a 
deeper bound state, and in the process enough energy 
is released such that both atom and dimer escape the 
trap. 

To connect the photoassociation rate to the mea¬ 
sured quantity we integrate the population evolution 
equation 


aV q 2 b,d q 2 b,t /.,.X 

where the dimer (trimer) photoassociation rate is 
multiplied by the numbers of pairs (trios) is the system. 
Since this experiment is clearly in the strong hierarchy 
regime, we use Eqs. (701 and 0 for the transition 
rates. 

To fix the value of L 2 we use Eq. ( 661 , where 
the relevant Feshbach resonance parameters are a^g « 
— 18.24ao and A « —237.8G [in]- Doing so, see Eqs. 
(721,(73), the transition rate is independent of A and 
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Table 1. The fitting parameters for the ^Li photoassociation 
experimental results. Here n is the number density and a is the 
broadening width. 


Case 

a(ao) 

n(10l2cm-3) 

Brf (G) 

Ea/h (MHz) 

cr (kHz) 

a 

806 

1.3 

0.6 

0.93 

25 

b 

629 

1.0 

0.52 

1.45 

20 


also of the effective charge (sq). 

The photon density is connected to the amplitude 
of the oscillating magnetic field i?rf through n^i = 
where is the vacuum permeability. 

Due to strong rf field there is power broadening 
which blurs the fine structure of the transition rate. 

To account for this effect we convolute the transition 
rates with a Gaussian, | 

O 

= -f=- [ duj'n^f{uj')e~^ 2 .^’ , (75) 

V/TTcr J 

where the Gaussian width a was fitted to the 
experimental results. 

The fitting parameters are shown in Table 

Our theory for photoassociation of ^Li ultracold 
atoms, as well as the experimental results of Ref. [S] are 
shown in Fig. To distinguish the power broadening 
effect the results without broadening are also shown 
in dotted curve. Our theory agrees well with the 
experimental results, and the fitting parameters are 
within the experimental uncertainty. However, the 
experimental resolution limits the ability to identify 
hne details of the association process. 



Figure 6. (Color online) The number of ^Li atoms remaining 
in the trap after photoassociation experiment, as function of the 
rf frequency. Shown are the experimental results from Ref. [^, 
for a/ao ~ 810 (upper plot) and a/ao ~ 627 (lower plot), as well 
as our fitting results. To distinguish the power broadening effect 
the results without broadening are also shown in dotted curve. 


7. Conclusion 

To conclude, we have applied the multipole expansion 
to study universal trimer photoassociation. One-body 
and two-body current were explored. For normal 
hierarchy case, where the one-body current dominates, 
the two leading s— and d— modes are studied and 
their relative contribution is shown to vary with 
temperature. In ultracold atoms, and other strong 
hierarchy systems, the relevant operator is found to 
be the two-body current. The log-periodic oscillations 
which modulate the cross section in both scenarios 
is found to be amplified for finite scattering length 
compared to the unitary point. We apply out theory 
to ^Li ultracold atoms photoassociation and show nice 
fit with the available experimental results. 
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Appendix A. Raynal-Revai coefficients 


Our calculation of the projector operator (34), provides 
an explicit formula for the = 0 or ly = 0 Raynal- 
Revai coefficients {Ixlyll'xly )kl ED) defined by 

= E(^Jy\Vy)KLylfj\^i^,)- (A.i) 

V I' 

I y 

Using the orthogonality of the hyperspherical harmon¬ 
ics we get, 

{lxly\l'Jy)KL = — --. (A.2) 

Putting lx = 0,ly = L and ai = 0, 

= (t)'^^(2z; + i)(2(; + i)x 


I'y L \ 

000) <^^^( 0 ) ■ 


(A.3) 
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Similar calculation for U = L,ly = 0 and ai = 7 r /2 
gives 

{L0\l'J'y)KL = {±fy^i2U, + l){2Uy + l)x 
I'y L \ 

0 0 0/ ‘/?fo(^/2) ' 

Appendix B. The hyperangular integrals 

The wave function is written as sum of Faddeev 
components, each in different coordinate set, 

1 ^ 

= -7='^YLo{yi)P!fia^)/sm2a,. (B.l) 

Our aim is to evaluate the hyperangular integral, 

In = J |^^r 2 ofe)cos 2 a, j $2 (^)(B. 2 ) 

Due to symmetry, one of the sums, say over j, can be 
dropped giving a factor of 3. The integral (B.2 1 can be 
further simplified utilizing the body-fixed coordinate 
system EO]. Consider the two body-fixed vectors Xio 
and yio such that iio • yio = 7 i- The transformation 
from the body-fixed to system to the laboratory system 
is given through the relation 


Ylmivi) — 'y ^ P ij.mi^)Yl^iyio), 


(B.3) 


where is the Wigner D-matrix, and the 

rotation w is the set of 3 Euler angles. Applying this 
transformation, for fixed a^. 


J dxi J dyiYijy^[yi)Yiijn'{yj) — 


.(B.4) 


J d'y^J2Yi*t,imo)Yi'y'iyjo) J dujD^*^{uj)D^y,^,{ui 

Now one can utilize the orthogonality of the Wigner 
D-matrix, 

J dajD^*rn{^)D^y,^,{uj) = (B.5) 

and the addition theorem 


07 _|_ 1 


(B. 6 ) 


where Pi(^ij) is the Legendre polynomials and 7 ij = 
ViQ ■ yjo- Choosing x^o = psm(at)z and y^o = 
pcos(ai)(^l - ilx P- ^iz), 


7y(a»,7i) = — 


COS ai 


± v^7 


'i COS ai 


\/l + 2 sin^ ai ± VSji sin 2ai 


(B.7) 


Finally, 


In = 


7^/2 




dji 


dai sin^ ai cos^ ai x 


<51 r 


i^k 


sin 2 Q!i sin 2ak 


where at = 0 ^( 01 , 71 ) is given by (14). 


(B. 8 ) 
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